home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / vopl / glvopl.lha / glvopl / src / spline.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-10-05  |  2.1 KB  |  117 lines
  1. #include "vopl.h"
  2.  
  3. /* 
  4.  *  spline
  5.  *
  6.  *     calculates a cubic spline through the data points
  7.  *     
  8.  */
  9. void 
  10. spline(x, y, np)
  11.     float     *x, *y;
  12.     int     np;
  13. {
  14.     float     *l, *u, *z, *a, *b, *c, *d, *h, *xx;
  15.     float    *alpha;
  16.     int     delta;
  17.     float     *newm1();
  18.     float     **newm2();
  19.  
  20.     float     t, xdivm, ax, ay, xdiv;
  21.     int     n, i, j;
  22.  
  23.     /* get memory for variables */
  24.  
  25.     h = newm1(np);
  26.     l = newm1(np);
  27.     u = newm1(np);
  28.     z = newm1(np);
  29.     a = newm1(np);
  30.     b = newm1(np);
  31.     c = newm1(np);
  32.     d = newm1(np);
  33.     xx = newm1(np);
  34.     alpha = newm1(np);
  35.  
  36.     n = np - 2;
  37.     for (i = 0; i <= n; i++) {
  38.         xx[i] = WhatX(x[i]);
  39.         h[i] = WhatX(x[i + 1]) - xx[i];
  40.         a[i] = WhatY(y[i]);
  41.     }
  42.  
  43.     n = np - 1;
  44.     xx[n] = WhatX(x[n]);
  45.     a[n] = WhatY(y[n]);
  46.  
  47.     if (plotdev.splinetype == CLAMPED) {
  48.         alpha[0] = 3 * (a[1] - a[0]) / h[0] - 3 * plotdev.s1;
  49.         alpha[n] = 3 * plotdev.sn - 3 * (a[n] - a[n - 1]) / h[n - 1];
  50.     }
  51.  
  52.     for (i = 1; i < n; i++) 
  53.         alpha[i] = 3 * (a[i + 1] * h[i - 1] - a[i] * (xx[i + 1] - xx[i - 1])  
  54.                + a[i - 1] * h[i]) / (h[i - 1] * h[i]);
  55.     
  56.     if (plotdev.splinetype == CLAMPED) {
  57.         l[0] = 2 * h[0];
  58.         u[0] = 0.5;
  59.         z[0] = alpha[0] / l[0];
  60.     } else {
  61.         l[0] = 1.0;
  62.         u[0] = 0.0;
  63.         z[0] = 0.0;
  64.     }
  65.  
  66.     for (i = 1; i <= n; i++) {
  67.         l[i] = 2 * (xx[i + 1] - xx[i - 1]) - h[i - 1] * u[i - 1];
  68.         u[i] = h[i] / l[i]; 
  69.         z[i] = (alpha[i] - (h[i - 1] * z[i - 1])) / l[i];
  70.     }
  71.  
  72.     if (plotdev.splinetype == CLAMPED) {
  73.         l[n] = h[n- 1] * (2.0 - u[n - 1]);
  74.         z[n] = (alpha[n] - h[n - 1] * z[n - 1]) / l[n];
  75.         c[n] = z[n];
  76.     } else {
  77.         l[n] = 1.0;
  78.         z[n] = 0.0;
  79.         c[n] = 0.0;
  80.     }
  81.  
  82.     for (j = n - 1; j >= 0; j--) {
  83.         c[j] = z[j] - (u[j] * c[j + 1]);
  84.         b[j] = (a[j + 1] - a[j]) / h[j] - (h[j] * (c[j + 1] + 2 * c[j]) / 3.0);
  85.         d[j] = (c[j + 1] - c[j]) / (3.0 * h[j]);
  86.     }
  87.  
  88.  
  89.     /*
  90.      * Draw the trace
  91.      */
  92.  
  93.     delta = (int)((float)plotdev.precision / (np - 1));
  94.     ay = a[0]; /*  + b[0] * t + c[0] * t * t + d[0] * t * t * t;*/
  95.     move2(xx[0], ay);
  96.  
  97.     for (j = 0; j <= n - 1; j++ ) {
  98.         xdiv = h[j] / (float)delta;
  99.         for (t = 0.0; t <= h[j]; t += xdiv) {
  100.             ax = xx[j] + t;
  101.             ay = a[j]  + b[j] * t + c[j] * t * t + d[j] * t * t * t;
  102.             draw2(ax, ay);
  103.         }
  104.     }
  105.  
  106.     free(h);
  107.     free(l);
  108.     free(u);
  109.     free(z);
  110.     free(a);
  111.     free(b);
  112.     free(c);
  113.     free(d);
  114.     free(xx);
  115.     free(alpha);
  116. }
  117.